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NOMENCLATURE 


a parameter with absolute value between 0 and 1 
coherence between ith and jth signals 

^ij^ii^jj 

spatial or time frequency (m"^ or sec'^) 
transfer function of filter between Xj and yj 
transfer function of filter between y^ and yj 
length scale of turbulence (m) 

Laplace transform s 

cross spectra between ith and jth signals 
friction velocity 

ith turbulent velocity component 

a parameter with absolute value between 0 and 1 

ith Gaussian white noise source 

Laplace transform of x ; 

lateral separation distance (m) 

ith simulated time series 

Laplace transform of y^ 

height above surface for ith signal (m) 

standard deviation of yj 



TECHNICAL PAPER 


ADDING COMPUTATIONALLY EFFICIENT REALISM TO MONTE CARLO 

TURBULENCE SIMULATION 

INTRODUCTION 


The goal of Monte Carlo turbulence simulation is to provide turbulence for flight simulation. 
The presence of turbulence in pilot-in-the-loop simulations increases the pilot’s workload as well as his 
sense of the realism of the simulation. A common complaint of pilots is that the turbulence provided 
is not realistic. The complaint is usually justified but not for the reason that the pilot believes. A real 
plane responds to turbulence in several ways. First, the random forces of lift and drag cause the center 
of gravity of the aircraft to translate in the three coordinate directions. Secondly, moments cause the 
craft to rotate about the roll, pitch, and yaw axes. Finally, the aircraft is an elastic body and the 
presence of aerodynamic forces and moments create a vibrational response in the airframe. The airframe 
“ringing” in many cases is a significant factor in creating the pilot perception of correct “feel.” The 
incorrect “feel” is unlikely to arise from inaccuracies in the stochastic properties of the turbulence 
which are provided but, rather from aspects of aircraft response which are not modeled. 

Of the three responses of the aircraft to turbulence mentioned above, only translation of the 
vehicle center of gravity is accounted for in common simulation practice. In some cases, only the vertical 
velocity component is simulated along the aircraft flight path. This turbulence component is typically 
simulated using explicit difference equations based on the Dryden spectral model, which is rational. 
Simulations comparing aircraft response to von Karman and Dryden turbulence were performed by Frost, 
et al. [ 1 ] and measurable differences were observed. These differences result from using less accurate 
rational models in generating one-dimensional turbulence. Roll and yaw moments cannot be realistically 
simulated unless variation of turbulent gusts across the body of the vehicle are simulated. Elastic airframe 
response cannot be simulated unless gust variation inputs are included. 

Several investigators have attempted to add three-dimensionality to turbulence simulation. Tatom, 
et al. [2] simulated turbulent gusts and gust gradients. This approach assumes that at any instant, turbu- 
lent velocities vary linearly in any direction across the body of the aircraft. The turbulence was generated 
using a non-dimensional von Karman model approach and simulation by this model involves manipulation 
of a sizable data base. The assumption of linear velocity variation is not realistic, even for an aircraft of 
moderate size, as was shown in the NASA Gust Gradient Program [3]. Even so, in the opinion of the 
author, this approach is much better than simple one-dimensional models. 

A highly accurate turbulence model which overcomes many conceptual difficulties was suggested 
by Campbell and Sandbom [4] . This technique combines remotely sensed atmospheric wind shear data 
with Monte Carlo generated three-dimensional turbulence. The problem with the model is that it requires 
the storage of a two-megabyte data base. Additionally, three-dimensional velocity cross correlations 
[R uv (x,y,z)] are not accounted for. At the present time, the technique would be useful only for large 

research simulators, but in the future as the cost of computer storage goes down, this or a derivative 
technique may be the method of choice. In any event, a long time will pass before any of the large air 
carriers will be willing to buy simulators with this kind of capability. 



The facts concerning turbulence simulation are that simulator manufacturers are willing to add 
any capability to their simulators that simulator users are willing to buy. The users cannot afford to 
purchase high fidelity turbulence simulation capability unless high computational efficiency and low 
storage requirements come with it. Even as computers become faster, cheaper, and more capable the 
situation will not change as fast as the computers. The reason is that flight simulation model software 
increases in size as better computers become available, and the requirement for more realistic turbulence 
must compete with other requirements. 

The purpose of this report is to document a concept which is an attempt to provide simulator 
users with three-dimensional turbulence simulation capability with a minimum increase in simulator 
computer demand. The concept is based on an earlier paper by Campbell [5], in which irrational von 
Karman spectra are approximated with rational functions. This approach is reviewed briefly in a later 
section, and its accuracy and computational advantages explained. The computational advantages of 
rational spectral models cannot be underestimated, and these will be explained in the next section. 


MONTE CARLO TURBULENCE SIMULATION 


The commonly applied engineering approach to Monte Carlo turbulence simulation is illustrated 
conceptually in Figure 1. Gaussian white noise is filtered by a linear filter whose filter function is chosen 
so that the spectrum of the output signal is the desired spectrum of turbulence. Because the filter is 
linear, the output signal is Gaussianly distributed. Some investigators use nonlinear filters (Reeves, et al. 
[6]) to better account for turbulent patchiness or other effects, but these techniques will not be 
considered here. 


GAUSSIAN 

WHITE 

NOISE 


FILTER 


"TURBULENCE” 


Figure 1. Monte Carlo turbulence simulation. 

Figure 1 is a conceptual diagram and many different approaches may be used to achieve the simu- 
lated turbulence. The filter function may take the form of von Karman, Kaimal, Dryden, or other 

spectra. The Kaimal and von Karman spectra are irrational having a falloff rate. The Dryden spectra 

fall off as F . If irrational spectra are used, a block of noise must be generated, transformed to the 
frequency domain, multiplied point by point by the irrational filter function, and then inverse trans- 
formed back to the time domain. Spectral models based on rational functions can be used to generate 
expheit finite difference equations, i.e. reccursion relationships. The latter approach is computationally 
much more efficient than the use of FFT’s for turbulence generation. For real time simulations where 
large blocks of mass storage are not available, the explicit difference equation approach is used almost 
exclusively. An understanding of the advantages of the rational versus the irrational approach is necessary 
to appreciate the technique proposed here, and the necessary explanations are given below. 

The equations governing incompressible flow are the Navier-Stokes equations and the equation 
of continuity. 
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( 1 ) 


1 '■y 

3uj/0t + Uj 3uj/3xj = - — 8p/3xj + v 3 z Uj/3xj3xj 
3uj/3xj = 0 (2) 


The Navier-Stokes equations are nonlinear partial differential equations for which very few analytical 
solutions are known. They can be solved numerically, but direct turbulent simulation of flows with the 
high Reynolds numbers of the atmosphere is not now feasible, nor will it be in the forseeable future. 
Given these facts, how can a realistic atmospheric simuland be created? The answer is to use one- 
dimensional measurements of turbulent spectra to generate differential systems which when random noise 
is input, not so random noise with the required spectra are output. Consider for a moment, the Dryden 
spectrum for the longitudinal velocity component (Table 1). The s-domain filter function corresponding 
to this spectrum is given by 

H L (s) = a 1 /(s + b 1 ) (3) 


Hj^(s) = Y(s)/X(s), where X(s) is the transform of the input white noise, and Y(s) is the transform of the 
filter output. Equation (3) can be rearranged as follows. 

(s + b } ) Y(s) = ajXfs) (4) 

If the inverse Laplace transform is performed on equation (4), the result is an ordinary, first-order 
differential equation with constant coefficients. 


dy/dt + bjy = ajx 


(5) 


Since bj is positive, this differential equation is stable since the filter function has only the one pole 
at -bj in the left half plane. Equation (5) can be discretized in many different ways, e.g. Runge-Kutta 

or predictor-corrector, but care must be taken that the resulting finite difference equations are also 
stable. Some commonly used discretizations do not satisfy this requirement, but Neuman and Foster 
[7] describe a scheme which guarantees that if the differential equation is stable the corresponding 
difference equation is stable. 

The longitudinal filter function for the von Karman spectra is given in Table 1 and if manipu- 
lated as above the result is 


(dj + s) 5 / 6 Y(s) = Cj X(s) 


( 6 ) 
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TABLE 1. SPECTRA, TRANSFER FUNCTIONS, AND DIFFERENTIAL EQUATIONS 
CORRESPONDING TO THE DRYDEN AND VON KARMAN MODELS 
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Before, Y was multiplied by s whereas in the above equation, it is multiplied by an irrational factor. 
Multiplying by s is equivalent to taking the derivative in the time domain. For the case of an irrational 
factor, calculation of derivatives of fractional order must be determined which is a cumbersome process 
[8], If derivatives of integer order are required, heuristically the irrational factor of Y(s) in equation 
(6) can be expanded in an infinite series. When transformed to the time domain, the result is a linear, 
ordinary differential equation of infinite order. Differential equations of infinite order cannot be 
discretized. 

The spectra, transfer functions, and differential equations (where possible) for the Dryden and 
von Karman spectra are tabulated in Table 1. For the rational Dryden spectra, first or second order, 
linear, nonhomogeneous differential equations with random input can be derived. The corresponding 
difference equations are of the following form [9] . 


y n +l = C 1 y n + c 2 y n -l + d l x n + d 2 x n-l 


(7) 


The calculation of an FFT requires a number of operations proportional to N k^N, where N 

is the number of points to be calculated. To generate N points of turbulence with an irrational spectrum 
requires the calculation of two FFT’s and a point-by-point multiplication by the filter function. The 
total number of operations is then proportional to 2N log^N + 4N multiplications. For the above differ- 
ence equation, the corresponding number of operations is 4N. If N = 1024, the number of multiplica- 
tions is six times greater for the FFT approach. For larger values of N, the computational advantage 
becomes correspondingly greater. The computational advantage is actually greater than this because for 
the above difference equation, only four coefficients must be calculated while considerably more sines 
and cosines must be calculated in the FFT evaluations. The disadvantage of the Dryden model is its 
inaccuracy. In fact, the Dryden model can be considered a crude approximation to the more accurate 
Von Karman model. Campbell [5] has developed a new method which applies accurate, rational approxi- 
mations to the von Karman spectra, and the result is simulated turbulence with rational spectra closely 
approximating the von Karman spectra. For the more complex transverse velocity component, 9N real 
multiplications are required to generate N points of turbulence. The computational efficiency is only 
9/4 as large as for the Dryden model, and still holds a strong advantage over the FFT approach. 

To this point, the only methods discussed that incorporate anything better than one-dimensional 
turbulence simulation were the studies by Tatom, et al. [2] and by Campbell and Sandborn [4] . The 
weaknesses of these two techniques were discussed. Other efforts which developed models for simulat- 
ing multiple time series should be mentioned. To match coherence, Fichtl, et al. [10] summed and 
phase shifted a series of random noise inputs. The coherence matching only becomes exact in the limit 
as the number of random inputs approach infinity. Hoshiya and Tieleman [11] developed a method for 
two time series. This method requires summations of large numbers of sine and cosine terms. Neither of 
the two models are as computationally efficient as desired. In the following section, a method is proposed 
which uses rational approximations, and multipoint turbulence simulation with cross correlation to 
achieve three-dimensionality. 

MULTI POINT TURBULENCE SIMULATION 


The problem addressed in this section is the generation of a number of time series with known 
spectra and cross^ spectra. Filter function forms will be developed for two series and for three series. 
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Generalizations to n dimensions will be apparent. The block diagrams will be developed in such a way 
that the form of the spectra and cross spectra of signals 1 and 2 will be the same whether only two 
signals or a hundred are being generated. Solution for the various filter functions in terms of the known 
spectra and cross spectra is thus simplified. Realizability of the required filter functions is another matter 
which cannot be presumed. For example, the product of a filter function and its complex conjugate must 
always be nonnegative. For some simulation situations negative values may be required to satisfy relation- 
ships between transfer functions and signal spectra and cross spectra. The question of realizability will be 
addressed briefly in this section. 

A block diagram which permits the generation of two time series with given spectra and cross 
spectrum is depicted in Figure 2 . The indicated filter functions are unknown, i.e. Hj, Hj2> and Ho. 

These are to be found in terms of Sjj, Sj2> and which are known quantities. Assuming that the 
Gaussian white noise inputs, Xj and X2> are independent and have unit variance, the three spectra are 
given by the following equations. 

Sn = IHjI 2 

S 12 = H 12 S 11 
S 22 -|S 12 I 2 /Sii+I h 2 I 2 

Equations (8), ( 9 ), and ( 10 ) when solved yield the required values for |Hjh, H^. and IH2 1 • Solutions 
for Hj and Hp are obvious. The solution for IH2 ! is more interesting. 

IH 2 I 2 - s 22 - is 12 l 2 /s n = S 22 (1 -coh 12 (0) • 01) 

'-y 

The question of realizability arises here in connection with |H->|~. Equation ( 11 ) reveals that a negative 
value of |H 2 r is not possible since S99 is positive and the coherence varies between 0 and 1. 



Figure 2 . Block diagram for generating two correlated time series. 


( 8 ) 

( 9 ) 

GO) 
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Figure 3 depicts the block diagram for the three time series case. In the figure all information 
flows from top to bottom, i.e., yj influences y 2 and y 3 , and ^2 influences y^, but yj is not affected 

by yg or y 2 and so on. The one way flow of information means that values of Hj, Hj 2 > and H 2 keep 

the same form that they had in the two signal case and as was given in equations ( 8 ) to (11). For the 
three signal case, only solutions of H j 3 , H 23 , and H 3 must be obtained. The results of the solutions 
are given below. 


H 13 = ( s i3 S 22 - s 12 S 23^ DEL (12) 

H 23 = ( S 11 S 23 - S 2 i s i 3 )/OEL (13) 

|H 3 I 2 = (S 1 1 S 22 S 33 + S 1 2 S 23 S 1 3 + S 13 S 21 S 32 “ S 11 > S 23 |2 “ S 33 |S 12 |2 " S 22 IS 13 I 2 )/DEL 

(14) 

where DEL = Sj JS 22 ~ 1 2 ^ 2 - DEL * s a l wa y s nonnegative. 



Figure 3. Block diagram for simulating three correlated time series. 

Realizability questions arise with regard to ]H 3 1 2 . The value of this quantity could conceivably 
become negative. The numerator of equation (14) is the determinant [S^]. Since equations ( 8 ) to (11) 
indicate that all the principle minors of rank one and two are nonnegative, a nonnegative value of the 
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determinant Sy means that the matrix would be positive definite. The matrix Sjj is Hermitian, i.e., 
Sjj = Sjj*. A nonrealizable example can be generated. 

1 x cx 

|H 3 I 2 = x 1 X / DEL (14) 

cx x 1 

The parameters, c and x have values between -1 and 1 . For a positive value of IH 3 I , the above deter- 
minant must be positive. The parameters are functions of frequency, but we are searching for negative 
values of the right side of equation (14) at only one frequency. The physical picture relating to equation 
(14) is that of three time histories recorded from an aircraft flying through isotropic turbulence (Fig. 4). 
The points are equally spaced so that Sp = S 23 = x. The cross spectra between the nonadjacent points, 

S 1 3 = cx. If the determinant in the above equation is evaluated, the result is given as follows. 

Det(Sjj) = 1 + 2x 2 c - 2x 2 - c 2 x 2 . (15) 



Figure 4. Physical example of simulation of three series 
at uniformly spaced points. 

For realizability, the determinant must be nonnegative. Regions of x-c space where the 
determinant is positive or negative may be specified by finding the curve where the determinant is zero. 
The above expression is quadratic in c and curves of zero value of the determinant are given by c = 1 /x 

and c = (2x“ - l)/x. Regions of negative value for equation (15) are indicated in Figure 5. Although the 
c-x plane is shown for large positive and negative values of x and c, in reality x and c only vary in 
absolute value between zero and one. The nonrealizable values in the c-x plane occur for x large, or for 
small separations of the three points of Figure 4 relative to the length scale of turbulence. Values of H 3 

were calculated for these small separations, and indeed computer errors related to negative square roots 
were encountered. Realizable transfer functions for larger separations are shown for different separations 
in the next section. 
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Figure 5. Three time series realizability in c-x space. 

The block diagrams shown in Figures 2 and 3 indicated a method of generating two and three 
time series with desired spectra and cross spectra. The extension to any number of series is straight 
forward and will not be presented here. Future research along these lines should be directed toward 
gaining an understanding of restrictions under which realizability is satisfied. Systems other than the 
ones presented here for generating correlated random signals might be discovered which assure 

realizability. 

In the preceding section, the Von Karman approximant (VKA) method [5] for approximating 
irrational spectra was described. In this section, forms of transfer functions suitable for generating 
correlated signals was described. The VKA method in combination with the transfer function forms 
presented here can provide a computationally efficient means for inserting three-dimensionality into 
Monte Carlo turbulence simulation. In the next section, example transfer functions for realistic situations 
will be shown. 
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EXAMPLE TRANSFER FUNCTIONS 


Turbulence encountered by an aircraft at altitude can be modeled as isotropic. Assuming the 
validity of Taylor’s frozen turbulence hypothesis and Von Karman turbulence, the cross spectra of the 
vertical velocity component was calculated by Houbolt and Sen [12]. Their expression for the cross 
spectrum is given in the following equation. 


(16) 


In the above expression, K a is the modified Bessel function of the second kind. Values of this function 

for a = 5/6 and 1 1/6 were obtained using a subroutine, BESSKA available from the National Center 
for Atmospheric Research (NCAR). Values for the transfer functions Hj, Hj2, H^, H^, H23, and H3 

are plotted in Figures 6 and 7. Figure 6 corresponds to a point spacing of y/L = 0.2, and 0.4. Figure 7 
is a corresponding plot for y/L = 0.05, and 0.1. For each of these separations H3 is realizable, but lor 

some smaller separations, it is not. These transfer functions are suitable for flight simulations at altitude 
where turbulence is more nearly isotropic. 
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Figure 6. Three series transfer functions for vertical velocities 
with spacing y/L = 0.2 and 0.4. 




Figure 7. Three series transfer functions for vertical velocities 
with y/L = 0.05 and 0.1. 

For flight at low altitude, an exponential form of the cross spectra is generally used. A nice 
summary of coherence functions in the lower atmosphere is given by Panofsky and Dutton [13]. The 
cross spectra used for demonstration purposes here is valid for the Cape Canaveral, Florida area. The form 
is taken from Fichtl, et al. [ 14] . 


S im = exp (-0.3465 Af/Af Q 5 + j27rgAf) (17) 

where Af - 27rf(z m /u m - z^). 

fg 5 = 0.00573 for the longitudinal component in strong winds. 

g = 0.7 for longitudinal component and the average of Zj and z m are less than 100 m. 

Zj = height above the ground. 
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For the purpose of demonstration, Zj = 10 m, Z2 = 20 m, and z^ = 30 m were selected. Using 
the logarithmic velocity profile and an assumed value of the friction velocity u* = 0.59, uj = 9.17 m/s, 
U2 = 10.18 m/s, and u^ = 10.79 m/s were calculated based on a typical roughness length of 0.02 m. 
Sjj was assumed to be in the form of the Von Karman spectrum for longitudinal velocity. Typical values 
of Lj = 21.4 m, L2 = 33.3 m, and L3 = 42.9 m were used. Gust standard deviations, a j = 1.185 m/s, 
02 = 1.517 m/s, and 03 = 1.75 m/s were calculated from information given by Tatom, et al. [2], The 

corresponding, realizable transfer functions are plotted in Figure 8. Interestingly, because of the large 
separations, and the fact that turbulent length scales increase as a fraction of height in lower planetary 
boundary layer (PBL), the transfer functions Hjj are very much smaller than the Hj. That being the case, 

for unrealizable solutions, the separations must be very small indeed to run into nonrealizability prob- 
lems. ’ 
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Figure 8. Three series transfer functions for longitudinal velocities and vertical separations 

(zj = 10 m, z-> = 20 nr, Z3 = 30 m). 




Another example which is interesting to consider is that of three dimensional, isotropic turbu- 
lence. The idea here is to generate three three-dimensional signals with given spectra and cross spectra. 
This is the most general simulation possible and would contain all one-dimensional cross spectra. The 
spectra in general must be in the following form. 


Sjj = E(f) (f 2 5y - f i fj)/47Tf 4 

(18) 

where f 2 = fj 2 + f 2 2 + f 3 2 and E(f) is a function which depends on the particular spectral model 
selected, i.e., Von Karman, Dryden, etc. Interestingly, the coherences cohy, and the transfer functions 
Hy are independent of the spectral model selected. The resulting transfer functions and coherences are 
given by the following equations. 

Hj = (E^irf 4 ) 1 / 2 (f 2 2 + f 3 2 )^ 2 

(19) 

H,2 * -flf 2 «f2 2 + f 3 2 > 

(20) 

H 2 = (E/47T1 4 ) 1 / 2 f f 3 /(f 2 2 + f 3 2 )^ 2 

(21) 

H 13 = - f l/ f 3 

(22) 

h 23 = _f 2/ f 3 

(23) 

h 3 = o 

(24) 

coh j 2 = (f! 2 f 2 2 )/( f 2 2 + f 3 2 ^ f l 2 + f 3 2) 

(25) 

cohis = (fj 2 f 3 2 )/(f 2 2 + f 3 2 )(f 1 2 + f 2 2 ) 

(26) 

coh 23 = (f 2 2f 3 2 )/( f l 2 + f 3 2 X f l 2 + f 2 2 ) 

(27) 


These equations indicate that the three-dimensional transfer functions are realizable for three- 
dimensional turbulence. In principle, it is possible to use these transfer functions in the system of Figure 
3. For isotropic turbulence, the resulting simulation would be the most realistic possible. 
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SUMMARY AND DIRECTIONS FOR FUTURE RESEARCH 


The purpose of this paper is to describe systems for generating multiple turbulence histories with 
given spectra and cross spectra. Using the systems depicted in Figures 2 and 3, two and three time 
histories can be generated. For the two time history case, the required transfer functions are always 
realizable. For the three time history case, small separations in space result in nonrealizable solutions. 
For realizable cases, once the exact transfer functions are determined, they may be fitted with rational 
functions as in the VKA method to generate sets of explicit difference equations which may be used to 
generate sets of turbulence time histories quickly and efficiently. Extensions to four or more time series 
is straight forward. 

Three examples of calculated transfer functions were presented. One was based on the Houbolt 
and Sen [12] cross spectral model for the vertical velocity component. For this model and small separa- 
tions, nonrealizable transfer functions occurred. For moderate separations, realizable transfer functions 
occurred. 

The second example model was for longitudinal velocities, separated vertically near the ground. 
Parameter values were selected as typical for the Cape Canaveral, Florida, area. Because of the strong 
diagonal dominance of the spectra tensor, no difficulties were encountered. 

The third example was for three-dimensional isotropic turbulence based on an arbitrary spectral 
model. This model yields realizable transfer functions for all situations. The cross spectra for this example 
were for different velocity components with zero separation. 

The system depicted in Figure 3 is not altogether satisfactory since it permits negative values of 
Hj [2]. The general problem of Monte Carlo simulation is to simulate a set of N time series with known 

spectra and cross spectra. If the outputs are to be Gaussian, then all filters are linear. The problem is to 
find the unknown filter functions from the known spectra and cross spectra. This concept is indicated 
in Figure 9. The block diagram of Figure 3 is a not quite successful attempt at realizing the system of 
Figure 9. Future research should be directed toward finding a realization of the general simulation 
problem. 

For the turbulence simulation problem, higher order realizable solutions to the general simulation 
problem should be found, and if irrational, fitted with stable rational functions as in the VKA method. 
The resulting stable rational functions can be used to generate stable, explicit difference equations for 
turbulence generation. These difference equations would be used to efficiently generate the turbulent 
velocities at several points over the three-dimensional body of an aircraft so that all aerodynamic forces, 
loads, and moments may be calculated, at least in principle. Only then will the atmospheric simuland be 
sufficiently accurate that flight simulation engineers can concentrate on other issues. 
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Figure 9. The general Monte Carlo simulation problem. 
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